Serratia (oligotyping)

Load packages, paths, functions

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot", "Biostrings", "openxlsx", "kableExtra")
invisible(lapply(packages, require, character.only = TRUE))

Preparation

Tables preparation

Seqtab

# move to oligotyping directory
setwd(paste0(path_oligo,"/serratia/oligotyping_Serratia_sequences-c1-s1-a0.0-A0-M10"))

# load the matrix count table
matrix_count <- read.table("MATRIX-COUNT.txt", header = TRUE) %>% t()

# arrange it
colnames(matrix_count) <- matrix_count[1,]
matrix_count <- matrix_count[-1,]
matrix_count <- matrix_count %>% as.data.frame()

# print it
matrix_count %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP41 S104 S25 S27 S30 S32 S36 S37 S38 S39 S44 S48 S49 S51 S77
A 10 24 10 7 442 3 10 3 4 8 0 26 2 2 0 4 3 2 1200 965 55 2132 46 666 17 2 43 13 314
G 4 6 1 0 207 0 7 2 2 5 1 5 2 1 40 0 5 0 655 298 36 874 7 306 6 0 23 5 0

Taxonomy

# move to oligotyping directory
setwd(paste0(path_oligo,"/serratia/oligotyping_Serratia_sequences-c1-s1-a0.0-A0-M10"))

# load the fasta table
fasta <- readDNAStringSet("OLIGO-REPRESENTATIVES.fasta")

# arrange it
fasta <- fasta %>% as.data.frame()
colnames(fasta) <- "seq"
fasta$oligotype <- rownames(fasta)
fasta <- fasta %>% dplyr::select(-c(seq))

# print it
fasta %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
oligotype
A A
G G

Change oligotype name by oligotype / MED nodes in the matrix count

# Reference file 

## move to tsv directory
setwd(path_tsv)

## load the reference table
ref_oligo_med2 <- read.table("2B_REF_info_serratia.tsv", sep="\t", header = TRUE)

## select only the 2 oligotypes of Serratia
ref_oligo_med2 <- ref_oligo_med2[!is.na(ref_oligo_med2$oligotype),]

## change order of columns
ref_oligo_med2 <- ref_oligo_med2 %>% select(c(seq, oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## create a column with reference name (will be used in plots)
ref_oligo_med2$ref <- paste0("oligotype_", ref_oligo_med2$OLIGO_oligotype_frequency_size, " / node_", ref_oligo_med2$MED_node_frequency_size)

## create a copy of fasta 
fasta2 <- fasta

# Matrix count

## create an oligotype column in the matrix count
matrix_count$oligotype <- rownames(matrix_count)

## change order of columns
matrix_count <- matrix_count %>% dplyr::select(c(oligotype, everything()))

## merge the matrix count and the reference dataframe
matrix_count2 <- matrix_count %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(c(oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size, ref, everything()))

## change rownames
rownames(matrix_count2) <- matrix_count2$ref

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(-c(oligotype, ref, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## print it
matrix_count2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP41 S104 S25 S27 S30 S32 S36 S37 S38 S39 S44 S48 S49 S51 S77
oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013 10 24 10 7 442 3 10 3 4 8 0 26 2 2 0 4 3 2 1200 965 55 2132 46 666 17 2 43 13 314
oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498 4 6 1 0 207 0 7 2 2 5 1 5 2 1 40 0 5 0 655 298 36 874 7 306 6 0 23 5 0 
## edit the fasta dataframe
fasta2 <- fasta2 %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")
rownames(fasta2) <- fasta2$ref
fasta2 <- fasta2 %>% dplyr::select(-c(MED_node_frequency_size, OLIGO_oligotype_frequency_size, oligotype))

## print it
fasta2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
ref
oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013 oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013
oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498 oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498

Metadata

metadata <- read.csv(paste0(path_metadata,"/metadata_22_06_21.csv"), sep=";")
rownames(metadata) <- metadata$Sample

Phyloseq object with oligotypes

# convert matrix_count into matrix and numeric
matrix_count <- matrix_count2 %>% as.matrix()
class(matrix_count) <- "numeric"

# phyloseq elements
OTU = otu_table(as.matrix(matrix_count), taxa_are_rows =TRUE)
TAX = tax_table(as.matrix(fasta2))
SAM = sample_data(metadata)

# phyloseq object
ps <- phyloseq(OTU, TAX, SAM)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 29 samples ]
## sample_data() Sample Data:       [ 29 samples by 13 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]
compute_read_counts(ps)
## [1] 8511
# remove blanks
ps <- subset_samples(ps, Strain!="Blank")
ps <- check_ps(ps)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 28 samples ]
## sample_data() Sample Data:       [ 28 samples by 13 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]

Create new metadata with Percent

Load ps with all samples (for final plot)

setwd(path_rdata)
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")
ps.filter <- check_ps(ps.filter)

Edit new metadata with Percent_serratia

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0)))

# add read depth in sample table of phyloseq object
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)

# select Wolbachia
ps.serratia <- ps.filter %>% subset_taxa(Genus=="Serratia")

# add read depth of Wolbachia
sample_data(ps.filter)$Read_serratia <- sample_sums(ps.serratia)
sample_data(ps.filter) %>% colnames()
##  [1] "Sample"         "Well"           "Strain"         "Field"         
##  [5] "Country"        "Organ"          "Species"        "Individual"    
##  [9] "Individuals"    "Date"           "Run"            "Control"       
## [13] "Dna"            "Species_italic" "Strain_italic"  "Read_depth"    
## [17] "is.neg"         "Read_serratia"
sample_data(ps.serratia) %>% colnames()
##  [1] "Sample"         "Well"           "Strain"         "Field"         
##  [5] "Country"        "Organ"          "Species"        "Individual"    
##  [9] "Individuals"    "Date"           "Run"            "Control"       
## [13] "Dna"            "Species_italic" "Strain_italic"  "Read_depth"    
## [17] "is.neg"
# add percent of Wolbachia
sample_data(ps.filter)$Percent_serratia <- sample_data(ps.filter)$Read_serratia / sample_data(ps.filter)$Read_depth

# round the percent of Wolbachia at 2 decimals
sample_data(ps.filter)$Percent_serratia <- sample_data(ps.filter)$Percent_serratia %>% round(2)

# extract metadata table
test <- data.frame(sample_data(ps.filter))

# merge this metadata table with the other
new.metadata <- data.frame(sample_data(ps)) %>% merge(test %>% dplyr::select(c(Sample, Read_depth, Read_serratia, Percent_serratia)), by="Sample")
new.metadata <- test[new.metadata$Sample %in% sample_names(ps),]
rownames(new.metadata) <- new.metadata$Sample

# print it
new.metadata %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
Sample Well Strain Field Country Organ Species Individual Individuals Date Run Control Dna Species_italic Strain_italic Read_depth is.neg Read_serratia Percent_serratia
CTC1 CTC1 G5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 29779 FALSE 14 0.00
CTC10 CTC10 D6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 77 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2609 FALSE 30 0.01
CTC11 CTC11 E6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 13874 FALSE 11 0.00
CTC12 CTC12 F6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 87 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1146 FALSE 7 0.01
CTC13 CTC13 G6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 18035 FALSE 649 0.04
CTC14 CTC14 H6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 101 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1708 FALSE 3 0.00
CTC15 CTC15 I6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 23180 FALSE 17 0.00
CTC2 CTC2 H5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 30692 FALSE 5 0.00
CTC3 CTC3 I5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 39920 FALSE 6 0.00
CTC4 CTC4 J5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 29 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2139 FALSE 13 0.01
CTC5 CTC5 K5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 15789 FALSE 1 0.00
CTC6 CTC6 L5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 19753 FALSE 31 0.00
CTC9 CTC9 C6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 67 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 4980 FALSE 3 0.00
NP14 NP14 K4 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti 1a 0 N/A run3 True sample 19 italic(“Aedes aegypti”) Field-Guadeloupe 7973 FALSE 0 0.00
NP2 NP2 K3 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 1c 0 N/A run3 True sample 89 italic(“Culex quinquefasciatus”) Field-Guadeloupe 648335 FALSE 0 0.00
NP20 NP20 E5 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti 3a 0 N/A run3 True sample 18 italic(“Aedes aegypti”) Field-Guadeloupe 136 FALSE 0 0.00
NP27 NP27 L5 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 7c 0 N/A run3 True sample 30 italic(“Culex quinquefasciatus”) Field-Guadeloupe 1234 FALSE 0 0.00
NP29 NP29 B6 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 9c 0 N/A run3 True sample 16 italic(“Culex quinquefasciatus”) Field-Guadeloupe 203 FALSE 0 0.00
NP30 NP30 C6 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 10c 0 N/A run3 True sample 24 italic(“Culex quinquefasciatus”) Field-Guadeloupe 228 FALSE 0 0.00
NP34 NP34 G6 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 14c 0 N/A run3 True sample 20 italic(“Culex quinquefasciatus”) Field-Guadeloupe 95 FALSE 0 0.00
NP35 NP35 H6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 7a 0 N/A run3 True sample 84 italic(“Aedes aegypti”) Field-Guadeloupe 196532 FALSE 0 0.00
NP36 NP36 I6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 8a 0 N/A run3 True sample 27 italic(“Aedes aegypti”) Field-Guadeloupe 249 FALSE 0 0.00
NP37 NP37 J6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 9a 0 N/A run3 True sample 52 italic(“Aedes aegypti”) Field-Guadeloupe 419340 FALSE 0 0.00
NP38 NP38 K6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 10a 0 N/A run3 True sample 62 italic(“Aedes aegypti”) Field-Guadeloupe 282479 FALSE 0 0.00
NP39 NP39 L6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 11a 0 N/A run3 True sample 47 italic(“Aedes aegypti”) Field-Guadeloupe 218684 FALSE 0 0.00
NP41 NP41 B7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 13a 0 N/A run3 True sample 93 italic(“Aedes aegypti”) Field-Guadeloupe 247152 FALSE 40 0.00
NP42 NP42 C7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 14a 0 N/A run3 True sample 85 italic(“Aedes aegypti”) Field-Guadeloupe 185157 FALSE 0 0.00
NP43 NP43 D7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 15a 0 N/A run3 True sample 104 italic(“Aedes aegypti”) Field-Guadeloupe 239335 FALSE 0 0.00
NP44 NP44 E7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 16a 0 N/A run3 True sample 50 italic(“Aedes aegypti”) Field-Guadeloupe 156879 FALSE 0 0.00
NP5 NP5 B4 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 2c 0 N/A run3 True sample 66 italic(“Culex quinquefasciatus”) Field-Guadeloupe 736159 FALSE 0 0.00
NP8 NP8 E4 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 3c 0 N/A run3 True sample 92 italic(“Culex quinquefasciatus”) Field-Guadeloupe 334799 FALSE 0 0.00
S100 S100 K7 Field - Camping Europe Field France Ovary Culex pipiens GL1 1 30/05/2017 run1 True sample 105 italic(“Culex pipiens”) Field-Camping~Europe 52486 FALSE 0 0.00
S102 S102 A8 Field - Camping Europe Field France Ovary Culex pipiens GL2 2 30/05/2017 run1 True sample 14 italic(“Culex pipiens”) Field-Camping~Europe 3456 FALSE 0 0.00
S104 S104 C8 Field - Camping Europe Field France Ovary Culex pipiens GL5 5 30/05/2017 run1 True sample 54 italic(“Culex pipiens”) Field-Camping~Europe 52403 FALSE 4 0.00
S105 S105 D8 Field - Camping Europe Field France Ovary Culex pipiens GL6 6 30/05/2017 run1 True sample 103 italic(“Culex pipiens”) Field-Camping~Europe 55577 FALSE 0 0.00
S106 S106 E8 Field - Camping Europe Field France Ovary Culex pipiens GL7 7 30/05/2017 run1 True sample 96 italic(“Culex pipiens”) Field-Camping~Europe 33053 FALSE 0 0.00
S107 S107 F8 Field - Camping Europe Field France Ovary Culex pipiens GL8 8 30/05/2017 run1 True sample 65 italic(“Culex pipiens”) Field-Camping~Europe 52154 FALSE 0 0.00
S108 S108 G8 Field - Camping Europe Field France Ovary Culex pipiens GL9 9 30/05/2017 run1 True sample 64 italic(“Culex pipiens”) Field-Camping~Europe 55735 FALSE 0 0.00
S109 S109 H8 Field - Camping Europe Field France Ovary Culex pipiens GL10 10 30/05/2017 run1 True sample 55 italic(“Culex pipiens”) Field-Camping~Europe 59023 FALSE 0 0.00
S110 S110 I8 Field - Camping Europe Field France Ovary Culex pipiens GL11 0 30/05/2017 run1 True sample 57 italic(“Culex pipiens”) Field-Camping~Europe 57377 FALSE 0 0.00
S121 S121 H1 Field - Bosc Field France Ovary Culex pipiens J26 22 28/06/2017 run2 True sample 37 italic(“Culex pipiens”) Field-Bosc 20361 FALSE 0 0.00
S122 S122 I1 Field - Bosc Field France Ovary Culex pipiens J27 23 28/06/2017 run2 True sample 52 italic(“Culex pipiens”) Field-Bosc 9803 FALSE 0 0.00
S123 S123 J1 Field - Bosc Field France Ovary Culex pipiens J28 24 28/06/2017 run2 True sample 102 italic(“Culex pipiens”) Field-Bosc 20130 FALSE 0 0.00
S124 S124 K1 Field - Bosc Field France Ovary Culex pipiens J29 25 28/06/2017 run2 True sample 69 italic(“Culex pipiens”) Field-Bosc 18146 FALSE 0 0.00
S126 S126 K6 Field - Bosc Field France Ovary Culex pipiens J30 26 28/06/2017 run2 True sample 99 italic(“Culex pipiens”) Field-Bosc 15235 FALSE 0 0.00
S127 S127 B2 Field - Bosc Field France Ovary Culex pipiens J31 27 28/06/2017 run2 True sample 38 italic(“Culex pipiens”) Field-Bosc 24696 FALSE 0 0.00
S128 S128 C2 Field - Bosc Field France Ovary Culex pipiens J32 28 28/06/2017 run2 True sample 71 italic(“Culex pipiens”) Field-Bosc 16305 FALSE 0 0.00
S146 S146 I3 Laboratory - Lavar Lab France Ovary Culex pipiens MW52 29 29/08/2017 run2 True sample 72 italic(“Culex pipiens”) Laboratory-Lavar 25012 FALSE 0 0.00
S147 S147 J3 Laboratory - Lavar Lab France Ovary Culex pipiens MW53 30 29/08/2017 run2 True sample 56 italic(“Culex pipiens”) Laboratory-Lavar 25171 FALSE 0 0.00
S148 S148 K3 Laboratory - Lavar Lab France Ovary Culex pipiens MW54 31 29/08/2017 run2 True sample 91 italic(“Culex pipiens”) Laboratory-Lavar 14164 FALSE 0 0.00
S150 S150 A4 Laboratory - Lavar Lab France Ovary Culex pipiens MW55 32 29/08/2017 run2 True sample 43 italic(“Culex pipiens”) Laboratory-Lavar 15081 FALSE 0 0.00
S151 S151 B4 Laboratory - Lavar Lab France Ovary Culex pipiens MW56 33 29/08/2017 run2 True sample 78 italic(“Culex pipiens”) Laboratory-Lavar 22944 FALSE 0 0.00
S152 S152 C4 Laboratory - Lavar Lab France Ovary Culex pipiens MW57 34 29/08/2017 run2 True sample 79 italic(“Culex pipiens”) Laboratory-Lavar 15082 FALSE 0 0.00
S153 S153 D4 Laboratory - Lavar Lab France Ovary Culex pipiens MW58 35 29/08/2017 run2 True sample 97 italic(“Culex pipiens”) Laboratory-Lavar 17040 FALSE 0 0.00
S154 S154 E4 Laboratory - Lavar Lab France Ovary Culex pipiens MW59 36 29/08/2017 run2 True sample 76 italic(“Culex pipiens”) Laboratory-Lavar 9626 FALSE 0 0.00
S160 S160 K4 Laboratory - Lavar Lab France Ovary Culex pipiens MW60 37 29/08/2017 run2 True sample 99 italic(“Culex pipiens”) Laboratory-Lavar 72508 FALSE 0 0.00
S162 S162 B5 Laboratory - Lavar Lab France Ovary Culex pipiens MW61 38 29/08/2017 run2 True sample 90 italic(“Culex pipiens”) Laboratory-Lavar 25180 FALSE 0 0.00
S163 S163 L6 Laboratory - Lavar Lab France Ovary Culex pipiens MW62 39 30/08/2017 run2 True sample 96 italic(“Culex pipiens”) Laboratory-Lavar 12333 FALSE 0 0.00
S164 S164 C5 Laboratory - Lavar Lab France Ovary Culex pipiens MW63 40 30/08/2017 run2 True sample 74 italic(“Culex pipiens”) Laboratory-Lavar 22368 FALSE 0 0.00
S165 S165 D5 Laboratory - Lavar Lab France Ovary Culex pipiens MW64 41 30/08/2017 run2 True sample 98 italic(“Culex pipiens”) Laboratory-Lavar 17731 FALSE 0 0.00
S166 S166 E5 Field - Camping Europe Field France Ovary Culex pipiens GL4 4 30/05/2017 run2 True sample 53 italic(“Culex pipiens”) Field-Camping~Europe 13979 FALSE 0 0.00
S167 S167 F5 Field - Bosc Field France Ovary Culex pipiens J32 28 28/06/2017 run2 True sample 59 italic(“Culex pipiens”) Field-Bosc 14048 FALSE 0 0.00
S169 S169 B7 Field - Camping Europe Field France Ovary Culex pipiens 5 43 16/05/2017 run2 True sample 95 italic(“Culex pipiens”) Field-Camping~Europe 11553 FALSE 0 0.00
S170 S170 C7 Field - Camping Europe Field France Ovary Culex pipiens 6 44 16/05/2017 run2 True sample 94 italic(“Culex pipiens”) Field-Camping~Europe 8852 FALSE 0 0.00
S18 S18 A1 Laboratory - Lavar Lab France Whole Culex pipiens MW75 0 30/08/2017 run1 True sample 3 italic(“Culex pipiens”) Laboratory-Lavar 4290 FALSE 0 0.00
S19 S19 B1 Laboratory - Lavar Lab France Whole Culex pipiens MW65 0 30/08/2017 run1 True sample 51 italic(“Culex pipiens”) Laboratory-Lavar 44527 FALSE 0 0.00
S20 S20 C1 Laboratory - Lavar Lab France Whole Culex pipiens MW66 0 30/08/2017 run1 True sample 39 italic(“Culex pipiens”) Laboratory-Lavar 42864 FALSE 0 0.00
S21 S21 D1 Laboratory - Lavar Lab France Whole Culex pipiens MW67 0 30/08/2017 run1 True sample 36 italic(“Culex pipiens”) Laboratory-Lavar 33798 FALSE 0 0.00
S22 S22 E1 Laboratory - Lavar Lab France Whole Culex pipiens MW68 0 30/08/2017 run1 True sample 54 italic(“Culex pipiens”) Laboratory-Lavar 19044 FALSE 0 0.00
S23 S23 F1 Laboratory - Lavar Lab France Whole Culex pipiens MW69 0 30/08/2017 run1 True sample 48 italic(“Culex pipiens”) Laboratory-Lavar 38172 FALSE 0 0.00
S24 S24 G1 Laboratory - Lavar Lab France Whole Culex pipiens MW70 0 30/08/2017 run1 True sample 70 italic(“Culex pipiens”) Laboratory-Lavar 42355 FALSE 0 0.00
S25 S25 H1 Laboratory - Lavar Lab France Whole Culex pipiens MW71 0 30/08/2017 run1 True sample 51 italic(“Culex pipiens”) Laboratory-Lavar 47688 FALSE 8 0.00
S26 S26 I1 Laboratory - Lavar Lab France Whole Culex pipiens MW72 0 30/08/2017 run1 True sample 17 italic(“Culex pipiens”) Laboratory-Lavar 5394 FALSE 0 0.00
S27 S27 J1 Laboratory - Lavar Lab France Whole Culex pipiens MW73 0 30/08/2017 run1 True sample 34 italic(“Culex pipiens”) Laboratory-Lavar 24558 FALSE 2 0.00
S28 S28 A2 Laboratory - Lavar Lab France Whole Culex pipiens MW74 0 30/08/2017 run1 True sample 6 italic(“Culex pipiens”) Laboratory-Lavar 4503 FALSE 0 0.00
S30 S30 K1 Laboratory - Lavar Lab France Whole Culex pipiens MW1 0 23/08/2017 run1 True sample 82 italic(“Culex pipiens”) Laboratory-Lavar 25353 FALSE 1855 0.07
S31 S31 L1 Laboratory - Lavar Lab France Whole Culex pipiens MW2 0 23/08/2017 run1 True sample 44 italic(“Culex pipiens”) Laboratory-Lavar 20417 FALSE 0 0.00
S32 S32 C2 Laboratory - Lavar Lab France Whole Culex pipiens MW3 0 23/08/2017 run1 True sample 22 italic(“Culex pipiens”) Laboratory-Lavar 12441 FALSE 1263 0.10
S33 S33 D2 Laboratory - Lavar Lab France Whole Culex pipiens MW4 0 23/08/2017 run1 True sample 26 italic(“Culex pipiens”) Laboratory-Lavar 33867 FALSE 0 0.00
S34 S34 E2 Laboratory - Lavar Lab France Whole Culex pipiens MW5 0 23/08/2017 run1 True sample 31 italic(“Culex pipiens”) Laboratory-Lavar 9367 FALSE 0 0.00
S35 S35 F2 Laboratory - Lavar Lab France Whole Culex pipiens MW6 0 23/08/2017 run1 True sample 23 italic(“Culex pipiens”) Laboratory-Lavar 11663 FALSE 0 0.00
S36 S36 G2 Laboratory - Lavar Lab France Whole Culex pipiens MW7 0 23/08/2017 run1 True sample 79 italic(“Culex pipiens”) Laboratory-Lavar 33020 FALSE 91 0.00
S37 S37 H2 Laboratory - Lavar Lab France Whole Culex pipiens MW8 0 23/08/2017 run1 True sample 88 italic(“Culex pipiens”) Laboratory-Lavar 18340 FALSE 3006 0.16
S38 S38 I2 Laboratory - Lavar Lab France Whole Culex pipiens MW9 0 23/08/2017 run1 True sample 63 italic(“Culex pipiens”) Laboratory-Lavar 54790 FALSE 53 0.00
S39 S39 J2 Laboratory - Lavar Lab France Whole Culex pipiens MW10 0 23/08/2017 run1 True sample 75 italic(“Culex pipiens”) Laboratory-Lavar 36273 FALSE 972 0.03
S40 S40 K2 Laboratory - Lavar Lab France Whole Culex pipiens MW11 0 23/08/2017 run1 True sample 83 italic(“Culex pipiens”) Laboratory-Lavar 44448 FALSE 0 0.00
S42 S42 A3 Field - Camping Europe Field France Whole Culex pipiens GLE1 0 30/05/2017 run1 True sample 1 italic(“Culex pipiens”) Field-Camping~Europe 4107 FALSE 0 0.00
S43 S43 B3 Field - Camping Europe Field France Whole Culex pipiens GLE2 0 30/05/2017 run1 True sample 11 italic(“Culex pipiens”) Field-Camping~Europe 9279 FALSE 0 0.00
S44 S44 C3 Field - Camping Europe Field France Whole Culex pipiens GLE3 0 30/05/2017 run1 True sample 4 italic(“Culex pipiens”) Field-Camping~Europe 8026 FALSE 23 0.00
S45 S45 D3 Field - Camping Europe Field France Whole Culex pipiens GLE4 0 30/05/2017 run1 True sample 13 italic(“Culex pipiens”) Field-Camping~Europe 18150 FALSE 0 0.00
S47 S47 F3 Field - Camping Europe Field France Whole Culex pipiens GLE6 0 30/05/2017 run1 True sample 15 italic(“Culex pipiens”) Field-Camping~Europe 1951 FALSE 0 0.00
S48 S48 G3 Field - Camping Europe Field France Whole Culex pipiens GLE7 0 30/05/2017 run1 True sample 60 italic(“Culex pipiens”) Field-Camping~Europe 56738 FALSE 2 0.00
S49 S49 H3 Field - Bosc Field France Whole Culex pipiens E1 0 28/06/2017 run1 True sample 28 italic(“Culex pipiens”) Field-Bosc 33498 FALSE 66 0.00
S50 S50 I3 Field - Bosc Field France Whole Culex pipiens E2 0 28/06/2017 run1 True sample 25 italic(“Culex pipiens”) Field-Bosc 28481 FALSE 0 0.00
S51 S51 J3 Field - Bosc Field France Whole Culex pipiens E3 0 28/06/2017 run1 True sample 40 italic(“Culex pipiens”) Field-Bosc 61788 FALSE 18 0.00
S52 S52 K3 Field - Bosc Field France Whole Culex pipiens E4 0 28/06/2017 run1 True sample 21 italic(“Culex pipiens”) Field-Bosc 21553 FALSE 0 0.00
S55 S55 B4 Field - Bosc Field France Whole Culex pipiens E6 0 28/06/2017 run1 True sample 46 italic(“Culex pipiens”) Field-Bosc 50447 FALSE 0 0.00
S56 S56 C4 Field - Bosc Field France Whole Culex pipiens E7 0 28/06/2017 run1 True sample 61 italic(“Culex pipiens”) Field-Bosc 42609 FALSE 0 0.00
S57 S57 D4 Field - Bosc Field France Whole Culex pipiens E8 0 28/06/2017 run1 True sample 86 italic(“Culex pipiens”) Field-Bosc 49157 FALSE 0 0.00
S58 S58 E4 Field - Bosc Field France Whole Culex pipiens E9 0 28/06/2017 run1 True sample 32 italic(“Culex pipiens”) Field-Bosc 30357 FALSE 0 0.00
S59 S59 F4 Field - Bosc Field France Whole Culex pipiens E10 0 28/06/2017 run1 True sample 33 italic(“Culex pipiens”) Field-Bosc 32798 FALSE 0 0.00
S60 S60 G4 Field - Bosc Field France Whole Culex pipiens E11 0 28/06/2017 run1 True sample 45 italic(“Culex pipiens”) Field-Bosc 44485 FALSE 0 0.00
S61 S61 H4 Field - Bosc Field France Whole Culex pipiens E12 0 28/06/2017 run1 True sample 41 italic(“Culex pipiens”) Field-Bosc 49545 FALSE 0 0.00
S63 S63 J4 Field - Bosc Field France Whole Culex pipiens E14 0 28/06/2017 run1 True sample 100 italic(“Culex pipiens”) Field-Bosc 53444 FALSE 0 0.00
S64 S64 K4 Field - Bosc Field France Whole Culex pipiens E15 0 28/06/2017 run1 True sample 81 italic(“Culex pipiens”) Field-Bosc 47628 FALSE 0 0.00
S79 S79 B6 Field - Camping Europe Field France Ovary Culex pipiens J16 12 28/06/2017 run1 True sample 68 italic(“Culex pipiens”) Field-Camping~Europe 59755 FALSE 0 0.00
S80 S80 C6 Field - Camping Europe Field France Ovary Culex pipiens J17 13 28/06/2017 run1 True sample 83 italic(“Culex pipiens”) Field-Camping~Europe 52788 FALSE 0 0.00
S83 S83 F6 Field - Camping Europe Field France Ovary Culex pipiens J20 16 28/06/2017 run1 True sample 73 italic(“Culex pipiens”) Field-Camping~Europe 42272 FALSE 0 0.00
S84 S84 G6 Field - Camping Europe Field France Ovary Culex pipiens J21 17 28/06/2017 run1 True sample 49 italic(“Culex pipiens”) Field-Camping~Europe 56676 FALSE 0 0.00
S85 S85 H6 Field - Camping Europe Field France Ovary Culex pipiens J22 18 28/06/2017 run1 True sample 35 italic(“Culex pipiens”) Field-Camping~Europe 41690 FALSE 0 0.00
S86 S86 I6 Field - Camping Europe Field France Ovary Culex pipiens J23 19 28/06/2017 run1 True sample 80 italic(“Culex pipiens”) Field-Camping~Europe 61984 FALSE 0 0.00
S87 S87 J6 Field - Bosc Field France Ovary Culex pipiens J24 20 28/06/2017 run1 True sample 58 italic(“Culex pipiens”) Field-Bosc 65958 FALSE 0 0.00
S88 S88 K6 Field - Bosc Field France Ovary Culex pipiens J25 21 28/06/2017 run1 True sample 106 italic(“Culex pipiens”) Field-Bosc 53102 FALSE 0 0.00 
# replace metadata in the created phyloseq object
sample_data(ps) <- sample_data(new.metadata)

Taxonomic structure

Count

col <- brewer.pal(7, "Pastel2")

# reshape data for plot
test3 <- test %>% select(c(Sample, Species, Strain, Organ, Read_depth, Read_serratia)) %>% reshape2::melt(id.vars=c("Sample", "Species", "Strain", "Organ"), vars=c("Read_depth", "Read_serratia"))

count_whole <- test3[test3$Organ=="Whole",]
count_ovary <- test3[test3$Organ=="Ovary",]

make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))

levels(count_whole$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_ovary$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_whole$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

levels(count_ovary$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))


# plot
p_count1 <- ggplot(count_whole, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=12, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=1.1), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
p_count2 <- ggplot(count_ovary, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
 facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
    ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=0.8), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
# afficher plot
p_count1
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_count2

# panels
p_group <- plot_grid(p_count1+theme(legend.position="none"), 
          p_count2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("B1", "B2"), c(0, 0), c(1, 0.5), size = 20)
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
legend_plot <- get_legend(p_count1 + theme(legend.position="bottom"))
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Whole (the most abundant nodes)

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0),
                                            nrow=2, byrow=TRUE))

# select whole
ps.filter.whole <- subset_samples(ps, Organ=="Whole")
ps.filter.whole <- prune_taxa(taxa_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole <- prune_samples(sample_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 26 samples ]
## sample_data() Sample Data:       [ 26 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]
# data pour plot
#data_for_plot2 <- taxo_data_fast(ps.filter.whole, method = "abundance")
data_for_plot2 <- taxo_data(ps.filter.whole, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot2$Name), "\",\n") %>% cat()
## "oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013.",
##  "oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498.",
data_for_plot2$Name <- data_for_plot2$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot2$Name <- as.factor(data_for_plot2$Name)
new_names <- c("oligotype_A (27) | size:6013 / N1114 (27) | size:6013.",
               "oligotype_G (23) | size:2498 / N1116 (23) | size:2498."
)

data_for_plot2$Name <- factor(data_for_plot2$Name, levels = new_names)

col_add <- brewer.pal(8, "Accent")

col <- c("oligotype_A (27) | size:6013 / N1114 (27) | size:6013."="#B136F5",
         "oligotype_G (23) | size:2498 / N1116 (23) | size:2498."="#B863FF")

levels(data_for_plot2$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot2$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

#data_for_plot2 <- data_for_plot2 %>% na.omit()

p2 <- ggplot(data_for_plot2, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p2

Ovary (the most abundant nodes)

# select ovary
ps.filter.ovary <- subset_samples(ps, Organ=="Ovary")
ps.filter.ovary <- prune_taxa(taxa_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary <- prune_samples(sample_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1 taxa and 1 samples ]
## sample_data() Sample Data:       [ 1 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 1 taxa by 1 taxonomic ranks ]
# data pour plot
data_for_plot3 <- taxo_data(ps.filter.ovary, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot3$Name), "\",\n") %>% cat()
## "oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013.",
data_for_plot3$Name <- data_for_plot3$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot3$Name <- as.factor(data_for_plot3$Name)

levels(data_for_plot3$Species)= c("Culex pipiens"=make.italic("Culex pipiens"))

levels(data_for_plot3$Strain) <- c("Camping~Europe")

p3 <- ggplot(data_for_plot3, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p3

# ps.filter.ovary@otu_table # only 1 oligotype
# 
# a <- data.frame(ps.filter.ovary@sam_data)
# b <- data.frame(ps.filter.ovary@tax_table)
# count_table <- data.frame(ps.filter.ovary@otu_table)
# 
# count_table %>%
#   kbl() %>%
#   kable_paper("hover", full_width = F)
# 
# a$Relative_Abundance <- 1.0000000000
# a$Name <- paste0(b$ref,".")
# 
# levels(a$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
#                "Culex pipiens"=make.italic("Culex pipiens"),
#                "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))
# 
# levels(a$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))
# 
# p3 <- ggplot(a, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+
#   geom_bar(position = "stack", stat = "identity")+
#   scale_fill_manual(values = "#B136F5")+
#   theme_bw() +
#   theme(axis.text.x = element_text(angle = 90)) +
#   ggtitle("") +
#   guide_italics+
#   theme(legend.title = element_text(size = 20), 
#         legend.position="bottom",
#         legend.text = element_text(size=14),
#         #legend.key.height = unit(1, 'cm'),
#         panel.spacing.y=unit(1, "lines"), 
#         panel.spacing.x=unit(0.8, "lines"),
#         panel.spacing=unit(0,"lines"),
#         strip.background=element_rect(color="grey30", fill="grey90"),
#         strip.text.x = element_text(size = 16),
#         panel.border=element_rect(color="grey90"),
#         axis.ticks.x=element_blank(),
#         axis.text.y = element_text(size=18)) +
#  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
#   labs(x="Sample", y="Relative abundance", fill="Oligotype / MED Node")
# 
# p3

Panels taxonomy of whole / ovary

legend_plot <- get_legend(p2 + theme(legend.position="bottom"))

# panels
p_group <- plot_grid(p2+theme(legend.position="none"), 
          p3+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A1", "A2"), c(0, 0), c(1, 0.5), size = 20)

p_taxo <- plot_grid(p_group, legend_plot, nrow=2, rel_heights = c(1, .1))
p_taxo

Save taxonomic plot

setwd(path_plot)

tiff("2Dg_OLIGO_counts_serratia.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
tiff("2Dg_OLIGO_taxonomic_serratia_whole.tiff", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("2Dg_OLIGO_taxonomic_serratia_ovary.tiff", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("2Dg_OLIGO_taxonomic_serratia.tiff", units="in", width=18, height=16, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("2Dg_OLIGO_counts_serratia_big.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("2Dg_OLIGO_counts_serratia_small.png", units="in", width=18, height=14, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("2Dg_OLIGO_taxonomic_serratia_whole.png", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("2Dg_OLIGO_taxonomic_serratia_ovary.png", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("2Dg_OLIGO_taxonomic_serratia_big.png", units="in", width=18, height=18, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("2Dg_OLIGO_taxonomic_serratia_small.png", units="in", width=18, height=14, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2

Make main plot

setwd(paste0(path_oligo,"/serratia/oligotyping_Serratia_sequences-c1-s1-a0.0-A0-M10/HTML-OUTPUT"))

img <- magick::image_read("entropy.png")
p_entropy <- magick::image_ggplot(img, interpolate = TRUE)
p_entropy+ theme(plot.margin = unit(c(-7,-2.5,-7,-0.5), "cm"))

p_entropy+ theme(plot.margin=unit(c(-7,-2,-12,-5), "mm"))

aligned <- plot_grid(p_taxo, 
                     p_counts, 
                     align="hv")

aligned

p_entropy2 <- plot_grid(p_entropy, nrow=1)+
  draw_plot_label(c("C"), c(0), c(1), size=20, hjust=-0.5)

p_entropy2

t_plot <- plot_grid(aligned, 
                    p_entropy2,
                    nrow=2, 
                    ncol=1, 
                    scale=1,
                    rel_heights=c(2,1))

t_plot

setwd(path_plot)

tiff("2Dg_OLIGO_main_serratia.tiff", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2
png("2Dg_OLIGO_main_serratia.png", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2